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ABSTRACT 

We have generalized a method for the numerical solution of hyperbolic systems of equations using 
a dynamic Voronoi tessellation of the computational domain. The Voronoi tessellation is used to 
generate moving computational meshes for the solution of multi-dimensional systems of conservation 
laws in finite-volume form. The mesh generating points are free to move with arbitrary velocity, with 
the choice of zero velocity resulting in an Eulerian formulation. Moving the points at the local fluid 
velocity makes the formulation effectively Lagrangian. We have written the TESS code to solve the 
equations of compressible hydrodynamics and magnetohydrodynamics for both relativistic and non- 
relativistic fluids on a dynamic Voronoi mesh. When run in Lagrangian mode, TESS is significantly 
less diffusive than fixed mesh codes and thus preserves contact discontinuities to high precision while 
also accurately capturing strong shock waves. TESS is written for Cartesian, spherical and cylindrical 
coordinates and is modular so that auxilliary physics solvers are readily integrated into the TESS 
framework and so that the TESS framework can be readily adapted to solve general systems of 
equations. We present results from a series of test problems to demonstrate the performance of TESS 
and to highlight some of the advantages of the dynamic tessellation method for solving challenging 
problems in astrophysical fluid dynamics. 

Subject headings: hydrodynamics - methods: numerical - relativity 



1. INTRODUCTION 

Many astrophysical gas dynamical systems involve the 
motion of compressible fluid with relativistic velocity or 
energy density. Techniques for the numerical solution of 
the equations governing the multi-dimensional dynamics 
of relativistic fluids have been developed greatly in recent 
years. Significant recent progress has been made in grid- 
based methods where the fuid is described using Eulerian 
meshes for hydrodynamics on both uniform meshes (see 
review by Marti & Muller, 2003, and referen ces therein) 



and with adaptive mesh refinement (AMR) ( Duncan fc 
Hughes|[T994l IZhang fc MacFadyen 1[200^ pander Hoist 



et al.||2008p, as well as using smoothed particle hydro- 



dynamics (Rosswog, 2010). Extensions of these methods 
have been implemented i ncluding for general relativistic 
magnetohydrodynamics |Font|200"3l |Gammie et al.|2003[ 
De Villiers fc Hawley||200^ |Del Zarma et al.|lj007| [Anni- 
nos et al. 120051 I Mignone et al.|2007| ITchekhovskoy et aE 
2007[ |Nagataki||2009[) and for dynamic spacetimes Q An- 
derson et al.||2006| |Cerda-Duran et al~||2008| |Etienne et 
al.l|2010p . 

in this paper we present a new method for relativistic 
gas dynamics based on a dynamic Voronoi tessellation 
of space. The Voronoi tessellation generates a numerical 
mesh which moves and distorts with an arbitrary velocity 
field depending on how the mesh-generating points in the 
tessellation are moved. While holding the mesh points 
fixed results in an Eulerian method, allowing them to 
move with the local fluid velocity results in an effectively 
Lagrangian method. In this paper we present the TESS 
code which we have written on a dynamic mesh gener- 
ated by successive Voronoi tesselations of the spatial do- 
main. TESS is a numerical code and general framework 
for both relativistic and non-relativistic implementations 
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of hydrodynamics and magnetohydrodynamics (MHD). 
The strength of the method is to retain the advantages of 
conservation-law formulation essential for accurate com- 
putation of relativistic gas dynamics, while gaining the 
advantages of a Lagrangian scheme. Of particular impor- 
tance, the mesh motion allows for contact discontinuities 
to propagate without the excessive numerical diffusion 
often found in fixed mesh computations. The preserva- 
tion of contact discontinuities is of great importance for 
problems involving the development of fluid instabilities 
and for reactive hydrodyamics where artifical diffusion 
of reactant species can lead to unphysical solutions. Us- 
ing mesh motion, TESS accurately solves the numericaly 
challenging case of relativistic shear flows (see Section 
3), a problem which underresolved Eulerian simulations 
calculate incorrectly (Zhang & MacFadyen, 2006). 

Lagrangian codes have had great success when em- 
ployed in one dimension, usually to treat problems with 
spherical symmetry. Multidimensional problems, how- 
ever, are more challenging for Lagrangian codes, due to 
distortion of the com putational me sh when complexities 
in the flow develop. Noh (1964) formulated a simple 



strategy for dealing with mesh distortion in multiple di 
mensions, which was to continuously remap the compu- 
tational domain as the system evolved, to prevent the 
mesh from becoming overly distorted. Codes employing 
this strategy are referred to as "Arbitrary Lagrangian- 
Eulerian" (ALE) codes. ALE codes solve the problem of 
mesh distortion, but a t the c ost of diffusive remaps. 



Borgers & Peskin 



(1987) addressed the problem of 
mesh distortion by using a tessellated computational do- 
main. The improvement was to employ the Voronoi tes- 
sellation; a unique decomposition of space into polyhc- 
dra based on a set of mesh-generating points (originally 
called "fluid markers"). The advantage of such an ap- 
proach is that the mesh-generating points can move freely 
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through the computational domain and can be added 
or removed to enable adaptive refinement. The Voronoi 
tessellation adjusts its shape and topology so that the 
c omputational mesh d oes not become overly distorted. 



Whitehurst (1995) also made use of domain tessel- 



lation tor mesh generation. His "signal method" was a 
finite- volume method, partially inspired by finite-element 
methods employed for solving problems with irregular 
boundaries. The method was first-order and conser- 
vative, and employed a Delaunay triangulation of the 
computational domain. The technique employed high- 
resolution shock-capturing techniques developed over the 
last decades for grid-based Eulerian codes, while also 
having the advantages that come from solving the fluid 
equations on a moving mesh. The "grid cells" in this case 
were Delaunay triangles, and the main source of diffusion 
came about during changes in the mesh topology. During 
such changes, the conserved quantities were redistributed 
evenly among the affected triangles. This process intro- 
duced diffusion similar to that encountered in ALE codes 
during remapping of the mesh, but because the diffusion 
was localized and only occurred during topology changes, 
the accumulated error was reduced. 



More recently, Springel (2010) developed a second 



order finite- volume approach to solving Euler's equations 
using a Voronoi tessellation of the computational do- 
main. This method was implemented in the AREPO 
code, which has recently been appli ed to simulate sta r 
formation in a cosmological context ( |Greif et al. ||201lj ). 
The advantage of using the Voronoi tessellation instead 
of the Delaunay triangulation is that the Voronoi cells do 
not suffer abrupt transformations during changes in mesh 
topology so that re-mapping and fluid re-distribution are 
uneccessary. This was the first time a second-order finite- 
volume Lagrangian method of this nature was proposed. 
Another important property of this method Galilean in- 
variance; the code's performance is independent of the 
velocity of the reference frame in which the calculation 
is performed. 

In this paper we generalize these developments to the 
case of relativistic hydrodynamics. We should not ex- 
pect to retain all of the advantages found in the Newto- 
nian case; in particular, the prospect of having a Lorentz- 
invariant formulation is not to be expected, since in stan- 
dard formulations mesh points are assumed to be simul- 
taneous at each timestep. 

The paper is structured as follows: In §2, we describe 
the formulation of the fluid equations, and the specific 
implementation in the code. In §3, we provide a series of 
test problems to determine the convergence rate and to 
illustrate what sort of problems naturally benefit from a 
Lagrangian approach. In §4, we demonstrate the code's 
usefulness in an astrophysical context, by looking at the 
Rayleigh- Taylor instability in a decelerating relativistic 
shock. In §5, we summarize our results. 

2. NUMERICAL METHOD 

We first give a brief description of the Voronoi tessella- 
tion, before explaining how it is used in solving the fluid 
equations. A tessellation is a decomposition of space into 
separate pieces, in our case polygons or polyhedra. For 
the moment, we restrict our attention to two-dimensional 
tesselations, mainly because they are easier to describe 
and to visualize. 



A Voronoi tessellation can be generated from some col- 
lection of points (see Fig. [T]). Each mesh-generating 
point (i.e. "mesh generator" ) will correspond to a poly- 
gon in the tessellation. The polygon associated with 
point P is defined to be the set of all points which are 
closer to P than to any other mesh generating point. 
Thus, if an edge is adjacent to the two Voronoi polygons 
of points P and Q, this edge lies in the line of points which 
is equidistant to P and Q. In general, we will refer to the 
polygons as "cells" and the edges as "faces" . This termi- 
nology makes it easier to generalize to arbitrary numbers 
of dimensions. Additionally, we will speak of the "vol- 
ume" of a cell, which in two dimensions will mean area 
and in one dimension will mean length. We will also re- 
fer to the "area" of a face, which in two dimensions will 
mean length, and in one dimension will just be a constant 
number which we set equal to unity. 

An important related tessellation is the Delaunay tri- 
angulation. Given a set of mesh generating points, one 
can generally form a triangulation with the mesh gener- 
ators as vertices. In a sense, the Delaunay tessellation is 
the "nicest possible" triangulation of the space, given by 
the following definition: for every triangle in the tessel- 
lation, the three vertices of the triangle all lie on some 
circle, C. In order for this to be a proper Delaunay trian- 
gle, no other mesh generators may lie within the circle C. 
This property almost uniquely defines the triangulation. 
In degenerate cases, where four or more mesh generators 
lie on a common circle, the triangulation is ambiguous for 
these mesh generators. Degenerate cases like this will not 
concern us greatly in this paper. 

For a given set of mesh generators, the Delaunay tessel- 
lation is the topological dual to the Voronoi tessellation. 
Every Delaunay edge corresponds to a Voronoi face, and 
every Voronoi vertex corresponds to a Delaunay triangle 
(in fact, the Voronoi vertex lies at the center of the circle 
generated by the three vertices of the Delaunay triangle). 
This fact will help us in constructing the Voronoi tessel- 
lation, because the problem can be reduced to finding 
the equivalent Delaunay triangulation. Since there is a 
straightforward test to check that a tessellation is Delau- 
nay, this is typically the easiest way to find out which 
mesh generators are neighbors. 

Before describing the details of how one might gener- 
ate a tessellation from a given set of points, it will be 
useful to write down the numerical formulation so that 
we know what information needs to be extracted from 
the tessellation procedure. We shall find that our for- 
mulation applies generally to any hyperbolic system of 
conservation laws, regardless of whether the underlying 
equations are relativistic. 

2.1. Finite Volume Formulation 
The equations being solved will always have the form: 

= s (i) 

U 



F = 



F 



In particular, for the case of Euler's equations: 
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U = | S l F J = SV + P8* 
E J \ Ev j + Pv j 



(2) 



(3) 
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Fig. 1. — The Voronoi (left) and Delaunay (middle) tessellations applied to the same set of mesh points. The right- most panel shows 
both tessellations superimposed, highlighting the duality between the two. The Voronoi tessellation is the one employed in TESS, but the 
Delaunay triangulation is used to determine neighbors in an intermediate stage of the tessellation process. 




Fig. 2. — The solid polyhedron traced out by a Voronoi cell 
during one timestep. The conservation law can be interpreted as 
"zero net flux" through all the faces of this spacetime polyhedron. 

where p is the density of the fluid, S is the momentum 
density, and E is the energy density. P is the pressure, 
and v is the flow velocity. For simplicity we consider the 
case where there are no source terms, S = 0. For the 
relativistic version: 

/ pu° \ ( pui \ 

U = phu l u° \F 3 = \ phu l u ] + PS ij . 

\ phu°u° - P- pu° J \ phu ^ - pv? J 

(4) 

Where now is the four-velocity. 

ph = p + e + P (5) 

and e is the internal energy density, which can be found 
from the equation of state: 

e = e(p,P). (6) 

For the case of an adiabatic equation of state, we have 

e = P/(T-l) (7) 

where T is the adiabatic index of the fluid. We consider 
general physical equations of state in the Appendix. To 
derive the finite volume form of these equations, we shall 
set the source term to zero for brevity, but the general- 
ization to a nonzero source term is straightforward. For 
concreteness, we shall work in 2+1 dimensions, but ev- 
erything said here easily generalizes to arbitrary num- 
bers of dimensions. Equation (flj) does not depend on 
the spacetime metric a priori, ana so we can write down 
a formulation independent of the metric. For the follow- 
ing derivation we shall assume a Euclidean metric, rather 
than a Minkowski metric. 



We now look at the evolution of a Voronoi cell over 
one timestep. It is assumed that the cell changes its 
shape and size by the linear motion of its faces, and that 
over a timestep it traces out a solid polyhedron in 2+1 
dimensions, whose top and bottom are surfaces of con- 
stant time (see Fig. (21). The shape is actually not quite a 
polyhedron, but this is an approximation we are making, 
which is valid in the limit of short timesteps. In practice, 
we resolve this issue and others by taking a high-order 
Runge-Kutta timestep. We integrate over this poly- 
hedral 2+1-dimensional volume P: 

J^rd 3 x = o. (8) 

We can easily convert this to an integral over the two- 
dimensional boundary of this solid: 

/ rn^d 2 x = (9) 

J dP 

where is the (euclidean) unit normal to the boundary. 
For the top and bottom of the polyhedron: 

n, = ( ^ ) • (10) 

The 2+1-dimensional unit normal to the other faces will 
be related to the 2-dimensional unit normal on a given 
timeslice, n, but it will also have a component in the time 
dimension because the face is moving with some velocity, 
w. If we assume that the face is not changing its size 
or orientation as it moves (another assumption which is 
resolved by taking a high-order Runge-Kutta timestep), 
it is straightforward to check that the 2+1-dimensional 
unit normal will be 

n„= 1 J-t™). (11) 

Now we evaluate the integrals in ^ by averaging the in- 
tegrated quantities over spacetime. In doing so, we need 
to know the 2+1 dimensional spacetime volume being in- 
tegrated over. For the top and bottom, this is straight- 
forward: 

f d 2 x = dV n , [ d 2 x = dV n+1 (12) 

J Bottom JTop 
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Fig. 3. — A space-time diagram of the Riemann problem with a moving face. The HLLC solver assumes two constant states in the 
interior of the Riemann fan, separated by a contact discontinuity. An eulerian code selects the state which houses the curve x/t = 0, 
whereas our formulation selects the state which houses x/t = w ■ h (in this example, the face lies in the region corresponding to U^,F^). 
Note that if a face moves superluminally, our scheme reduces to upwinding. 



where dV is the cell volume at the beginning or end of 
the timestep. For the other faces, it is easy to check that: 



d 2 x = y/l + (w- n) 2 dAAt. 



(13) 



Face 



Recall that "dA" is the face "area", so it refers to the 
length of a Voronoi edge. Note that our factors of 
\J\ + (w ■ h) 2 will end up cancelling, which is to be ex- 
pected if our formulation is independent of the spacetime 
metric. If we interpret U n as the cell-averaged value of U 
at timestep n, and let Fij denote the time-averaged flux 
through the face adjacent to cells i and j, and likewise 
use Uij to denote the time-and-area-averaged value of U 
on this same face, we get the following result: 



= J dP f> 1 n fi d 2 x = 



(14) 



jjn+ljyn+l _ u n dV n + At £ (F tJ - w ■ hU l3 )dAj 

cellj 

This gives us a simple prescription for how to evolve the 
conserved variables from one timestep to the next, as- 
suming we know the time-averaged fluxes and conserved 
quantities on the faces, and the velocity of each face: 



U n+1 dV n+1 = U n dV n 



cellj 



Wij ■ hUij)dAj 



(15) 

Of course, this result is not surprising at all. The pre- 
scription merely tells us to add an advective term to the 
flux (F — > F — wU), and evolve things in a way analo- 
gous to the fixed-grid approach. It is worth noting that 
our formulation does not depend on the physical content 
of the equations expressed by ([I]). In particular, it does 
not matter whether the velocity w exceeds the speed of 
light. In order for this to be possible, we must be careful 
about our physical interpretation for the mesh itself. For 
example, it is not necessarily meaningful to speak of the 
"rest frame" of a cell or face. 

2.2. Riemann Solver 



Equation (15) approaches an exact result, in the limit 



of short timesteps. This means that all of the numerical 
error due to the spatial discretization must stem from our 
estimates of the time-averaged values of the fluxes and 



conserved quantities on the faces of the Voronoi cells. In 
TESS, we estimate these fluxes using an approximate 
Riemann solver. This means we assume that, at the 
beginning of the timestep, there is a uniform state on 
either side of the face. Generally speaking, a Riemann 
solver estimates the time-averaged flux through the face 
by evolving this constant-state initial condition, cither 
exactly or approximately. 

It turns out that getting the most out of the La- 
grangian nature of our scheme is highly dependent on 
how we solve the Riemann pr oblem on each ce ll face. The 
Riemann solver in AREPO ( [SpringeF 2010[ ) is effective 
because it is Galilean-invariant; the Kiemann problem is 
solved in the rest frame of the face. In our case we cannot 
assume it makes sense to even speak of a "rest frame" for 
a face, so we cannot have a method which is Galilean- 
invariant (this is to be expected anyway for relativistic 
problems). Nonetheless, we can still have a code which 
preserves contact discontinuities to very high accuracy. 
What we require is a Riemann solver which is exact for 
a pure advection problem. The HLLC Riemann solver, 
appropriately employed, meets this requirement. In the 
TESS code, we employ t he HLLC solver for relati vistic 
problems as described by Mignone & Bodo ( 2005 ) . 

A pure advection problem lias constant pressure and 
velocity on both states, and two different densities. The 
HLLC solver divides spacetime into four pieces: the orig- 
inal left and right state, and two interior states known 
as *L and *R, separated by a contact discontinuity. For 
advection, we have the following: 



P*l = PL 
P*R = PR 
P*L = P*R = P L =P R 

v*l = v* R = vl = vr. 



(16) 



This is the exact solution to the advection problem, and 
hence when using the HLLC solver, we can solve the ad- 
vective Riemann problem to machine precision, as long as 
the correct starred state is chosen, i.e. the solver should 
choose the spacetime region which houses the path traced 
out by the face as it moves with velocity w (see Fig. [3]). 
The values F* and U* found in this spacetime region 
are t he values used in the update step given by equa- 
tion (15). The reason that the HLLC solver is so crucial 
is that it accurately calculates advective fluxes. Since 
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the numerical error caused by the spatial discretization 
is entirely housed in estimating the time-averaged fluxes 
F ■ n — w ■ nU, and if the velocity of each face is very 
close to the fluid velocity, then for HLLC the advective 
fluxes will cancel, so that the numerical error for advec- 
tive fluxes will be small. In the particular case of pure 
advection, advective fluxes completely cancel, meaning 
that the advection problem is solved exactly to machine 
precision. This property is extremely important for pre- 
serving contact discontinuities. 

2.3. Primitive Variable Solver 

The Riemann solver takes in two states correspond- 
ing to two adjacent Voronoi cells. To use the informa- 
tion in a cell to solve the Riemann problem for a given 
face, we need the primitive variables (density, pressure, 
and four velocity) on either side of the face. In order to 
find these primitive variables, we need to invert formu- 
las ([3]) or Q for the conserved variables on either side. 
This can be done using a Newton-Raphson rootfinding 
scheme. For relativistic hydrodynamics with an adia- 
batic equation of state, this solver is not difficult to 
write. For an arbitrary equation of state, we incorporate 
the temperature as an additional variable which must 
be solved for. In addition to an arbitrary equation of 
state, we want TESS to become a very general code ca- 
pable of solving wide classes of problems, for example the 
equations of general relativistic magnetohydrodynamics 
(GRMHD). This makes the Newton-Raphson step some- 
what more complicated, but this has already been em- 
ployed in TESS even though we have not yet employed 
hydrodynamic variables which would necessitate such a 
solver. The GRMHD prim itive variable solv er is based 
on the method described by |Noble et al. |( [2006[ ) and uses a 
three-dimensional Newton-Raphson step, solving for the 
variables W = j 2 ph, K — j 2 — 1, and the temperature, 
T (see the Appendix for more details). 

2.4. Reconstruction of Face- Centered Primitive 
Variables 

The root finding method described above determines 
the primitive variables (density, pressure, and four- 
velocity) at the cell centers. To solve the Riemann prob- 
lem on each face, we must extrapolate the primitive vari- 
ables from the cell centers to the face centers. If we 
assume all our variables to be piecewise constant, then 
we can assume they have the same value on the face as 
they do at the center of mass of a cell. However, if we 
want accuracy higher than first order in space, we need 
to extrapolate the variable values based on the values 
in neighboring cells. We use piecew ise linear reconst ruc- 
tion, using the method derived by Springel (20101 for 
calculating the variable gradients in each cell. We repeat 
the results here. Assume we have some variable, 0, for 
which we would like to calculate the gradient at cell i 
based on the values at adjacent cells. The formula we 
use is 



I % V i I in Z / in 



)• (17) 



We use this gradient to extrapolate primitive variable 
values via: 



0(/tf) = 0i + (v0}.-(4--s i ). 



Here, r represents the location of a mesh generator, s 
represents the location of the center of mass of a cell, 
and ftj is the center of mass of the face adjacent to cells 
i and j . Note that primitive variables are defined at s, not 
at the mesh generators. This prescription would lead to a 
code which is second order in space, but it is well known 
that piecewise linear reconstruction can cause numerical 
oscialltions in the neighborhood of shocks. To deal with 
this, we need to constrain the estimated gradients in the 
neighborhood of sharp discontinuities. In other words, 
we need to construct a generalization of the slope limitcrs 
used in grid-based Eulerian codes. 

AREPO uses a slope limit er which could be considered 
a generalization of minmod ( Kurganov fc Tadmor |2"000| , 
but one which does not have the total variation dimin- 
ishing (TVD) property. As a result, this slope limiter 
caused mild oscillations in some calculations. TVD is an 
especially important property for relativistic hydrody- 
namics, since oscillatory behavior in the conserved vari- 
ables can potentially cause wild variation of the primitive 
variables, particularly in situations with large Lorentz 
factors. To address this problem, we optionally employ 
an alternative slope limiter which is much more conser- 
vative. AREPO uses the slope-limited gradient, 



= a, ( V(j> 

i \ 

min(l, tpij) 



ipi 



Ai/Jij < 
Aipii = 



(19) 
(20) 



(21) 



Our method is similar, but replaces (21 1 with 



( max(6('il'j 
tpij — 1 max(9(ipj 
1 



^)/A^,0) 
^)/A^-,0) 



N>a > 

Ai/iij < (22) 
A^i,- = 



where 9 is generally set to unity, but reduced if a still 
more diffusive scheme is desired. The slope limiter in 



( 22 ) enfo rces monotonicity, though it is more diffusive 
than ( |21[ ). In practice, we find it to be much more robust, 
so it has typically been employed in problems involving 
strong shocks. 

2.5. Time Integration 

Our time integration is based on the method of lines, 
performing a Runge-Kutta timestep for the time evolu- 
tion of the conserved variables, and for the motion of 
the mesh points. For most problems, we use a third or- 
der Runge-Kutta timestep which is TVD (total variation 
diminishing,) and which updates both the values of the 
conserved variables, and the positions of the mesh gen- 
erating points. The timestep is Courant-limited: 



At = C c ji ■ min{ 



Ri 



-) 



(23) 



(18) 



C c fi is the Courant factor, typically chosen between 0.2 

and 0.5. Ri is the effective radius of a cell, R = dV/n 
in 2D. A™ oa: is the eigenvalue in cell i with the largest 
magnitude. Currently the code is set up to operate with 
a single global timestep. To take a timestep from U n to 
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U n+1 which is third-order in time, for example, we use 
the following prescription: 

{/(!) = U n + AtL([/",r tt ) 
fW = f^+Atw 71 



C/(2) = |C/« + I[/(D + iAiL(C/( 1 ),r< 1 )) 



f(2) = | p» 



4*1 
4' 



. 2fj(2) 

2f(2) , 
3 T 



(24) 



|-|AiL([/( 2 ),r< 2 )) 



Here, L is an operator representing the numerically 
integrated time derivative of U, and the variables 
[A 2 ) } and f 1 - 2 ' represent intermediate states in 
the time integration. 

2.6. The Voronoi Mesh 



Equation (151 tells us exactly what geometric informa- 



tion we need m order to evolve the conserved quantities. 
We need to know the following about the Voronoi cells: 

• which cells are neighbors 

• the volume of each cell 

• the area of each face 

• the velocity of each face 

• the center of mass of each cell 

• the center of mass of each face 

The last two elements of this list are necessary for the 
piecewise linear reconstruction of primitive variables. We 
must determine how to extract all of this information 
given the positions and velocities of all the mesh gen- 
erating points. The velocities of mesh generators are 
freely specifiable, and we shall typically choose to set 
them equal to the local fluid velocity. 

Before we determine this completely, we can show that 
all of the above can be calculated easily if we know which 
cells are neighbors, and if we also know the center of mass 
and area of each face. In other words, when performing 
the Voronoi tessellation, this will be the only information 
we need to store. 

Given a single mesh generator and its neighbors, it is 
straightforward to calculate its cell's volume, if the area 
of each face is known. This is done by dividing the cell 
into pyramids (in 2D, triangles), each of whose tip is the 
mesh generating point, and whose base is a Voronoi face. 
Then the volume of the cell is the sum of the volumes of 
all the pyramids: 



dVi = EjV Aj 



(25) 



The volume of a pyramid can be expressed generally in 
D-dimcnsions: 



Va = (area of base) (height )/D 



(26) 



Because the face is in the plane halfway between the two 
mesh generating points, the height should be half the 
distance between these points. 



V Aj = dA ij (\r ij \/2)/D 



(27) 



D 



(28) 



Similarly, the center of mass of a cell can be directly 
calculated from the area and center of mass of the faces. 
We can use the weighted average over pyramids: 



Si — -jyT,jVAjSAj 



1 ^ dAj(\ rij \/2) 



dV, 



D 



-S A j- 



(29) 
(30) 



The center of mass of a pyramid also depends on the 
number of dimensions: 



S A j 



D 



D + l 



fij + 



1 



D + l 



(31) 



Here, fij denotes the center of mass of the face adjacent 
to cells i and j. 

Next, we need to determine the velocity of the faces. It 
is assumed here that the mesh-generating points them- 
selves have been given some velocity Wi, typically the 
local fluid velocity (though corrections can be added 
to steer t he cells in ways that make the mesh better- 



behaved). Springel | (120101) showed that the velocity of 
a face can be calculated from the position and velocity 
of the mesh generating points and the center of mass of 
the face. The result is the average of the velocity of the 
two adjacent mesh generators, plus a correction due to 
the fact that the center of mass of the face is not gen- 
erally at the midpoint between the two mesh generating 
points, and so acquires a velocity due to rotation about 
this midpoint: 



w F — (w L + w R )/2 + w 



?' = (w L - w R ) ■ (/ - (r L + r R )/2) 



r R - r L 



(32) 
(33) 



Again, f is the center of mass of a face. We can use equa- 



tions ( 25 - 33 ) to pare down the list of information that 



we neea to extract directly from the tessellation itself. 
The tessellation procedure now only consists in deter- 
mining the following: 

• which cells are neighbors 

• the area of each face 

• the center of mass of each face 

All relevant geometrical information can be easily ex- 
tracted from this data. This is advantageous because 
it means the tessellation can take up a relatively small 
amount of memory. The tessellation procedure consists 
of generating a new set of faces each timestep, based on 
the locations of the mesh generators. 

In one dimension, the tessellation procedure is trivial; 
neighboring cells do not change, the face area is always 
set to unity, and the face center of mass is simply the 
midpoint between the two mesh generators. The two- 
dimensional version turns out to be surprisingly sim- 
ple, because we use the tessellation from the previous 
timestep to generate the new faces. Since we know the 
neighbors of each cell on the previous timestep, we can 
use the neighbors of neighbors ("friends of friends") of 
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Fig. 4. — A single step of the tesselation process. We decide 
whether to keep point B by asking whether point C is outside of 
the circle generated by X, A, and B. 



a cell as a list of candidates for the neighbors on the 
next timestep. Because the length of each step in time 
is Courant-limited, the tessellation will not change sig- 
nificantly in one timestep, and hence this list of candi- 
dates is big enough for our purposes. Optionally, we 
can choose to use "neighbors of neighbors of neighbors" 
but we have not found this to make a difference in any 
scenario we've encountered. Already having this list of 
candidates simplifies the algorithm immensely, as in prin- 
ciple it can reduce to a brute force guess-and-check pro- 
cedure using this small list of candidates. In practice, the 
2D algorithm is not totally brute force, but very simple 
nonetheless. 

We consider a single mesh generating point X and all of 
its "friends of friends" . First we find the nearest neighbor 
to X (which is guaranteed to share a face with X). Call 
this neighbor Y. Next we take the rest of potential neigh- 
bors and order the list by clockwise angle with respect to 
the segment XY. What follows is an elimination process; 
at the end of this process, the elements in the list will 
be exactly those which share a face with X. At each step 
in the process, we consider three consecutive elements of 
the list; call them A,B,C (see Fig. We denote the 
element before A "P", and the element after C "Q". It 
is determined whether or not to keep point B in the list, 
by checking whether C lies within the circle generated 
by X, A, and B. If it does, then we remove B from the 
list, and take one step backward (checking whether C lies 
within XPA). If C does not lie within the original circle, 
we keep B and move forward to check triangle BCQ). 
More concisely, if Qxab contains C, remove B from the 
list and take one step back: A — > P. B — > A. Otherwise 
take one step forward: C->Q,B->C,A^B. Fig. [5] 
demonstrates an example of this full procedure. 

Note that this algorithm is not sensitive to the presence 
of degenerate sets of points (that is, sets of four points 
which all lie on the same circle). For practical purposes, 
it will not matter whether our code chooses to accept or 
reject a point in this configuration, because a degeneracy 
in the Delaunay tessellation corresponds to a face of zero 
area in the Voronoi tessellation. If the face has zero area, 
then there will be zero flux through it, and hence it will 
have no influence on the resulting physics. 



Once this operation has been performed for all mem- 
bers of the list (i.e., point B is now point Y, the near- 
est neighbor) all remaining list members are neighbors 
of point X. All that remains is to calculate the areas 
and centers of mass of the faces, which is straightforward 
given the vertices of the polygon generated by this list. 
These vertices are the centers of the circles generated by 
consecutive triples in the list. 

The details of the tessellation algorithm do not com- 
pletely extend to three dimensions, but the idea of using 
friends-of-friends as a list of candidates still works, so in 
the worst case scenario we could use a brute force guess- 
and-check algorithm when D=3. One might ask whether 
there is a major efficiency advantage to this tessellation 
procedure over more conventional ones such as direct in- 
sertion. While it is not clear which method should be 
the fastest (given an optimized implementation) , it is as- 
sumed they are comparably efficient. Moreover, while 
the tessellation procedure takes up a non-negligible per- 
centage of the code's overall runtime, it does not take up 
a majority of the runtime, and as such there is no major 
incentive to optimize its efficiency. The main advantage 
to the method described here is that the algorithm is 
very simple and does not require making a lot of ex- 
ceptions. Additionally this algorithm is expected to be 
very easy to parallelize, because the tessellation is per- 
formed locally. In the most simple prescription, we could 
make the code parallel via a simple domain decomposi- 
tion, where different processes only share boundary data. 
The main disadvantage to the tessellation procedure is 
that we must have an approximate tessellation to begin 
with, so that we can use "friends of friends" for our ini- 
tial pool of neighbors. In principle, this simply amounts 
to saving the tessellation from the previous timestep, but 
in practice this makes the implementation of boundary 
conditions a bit more complicated. Not only do we need 
to create an appropriate set of "ghost cells" outside the 
boundary, but we need to generate "ghost faces" as well 
(the tessellation procedure won't automatically do this 
for us, because it relies on having an approximate tessel- 
lation from a previous timestep). This also adds some 
complication to the parallelization of the code, for the 
same reasons. 

One might wonder why we worry so much about the 
positions of ghost cells. For periodic boundary con- 
ditions, we don't really need ghost cells at all, since 
we could in principle associate left-most neighbors with 
right-most neighbors and so on, so that no boundary 
need be created. For reflecting boundaries, one might 
hope we could have a fixed set of ghost cells lining our 
reflecting wall, so that we don't need to generate a new 
set each timestep. Unfortunately, we only have control 
over the positions of mesh generators; we don't directly 
control shape of the cells. Thus, if we want a flat wall, 
we need to have the mesh generators reflected across the 
wall. The only reason we include ghost cells for periodic 
boundary conditions is for the sake of overall simplic- 
ity in the code (fewer parts of the code depend on the 
choice of boundary conditions). In this case, ghost cells 
are translated to the opposite end of the domain with 
respect to their "real" counterparts. 

The boundary conditions are set in each dimension se- 
quentially. The first step involves flagging cells according 
to whether they are inside or outside the boundary. If we 
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Fig. 5. — A demonstration of the tessellation algorithm in two dimensions. Potential neighbors are sequentially tested and either accepted 
or rejected based on the criterion in Fig. [4] Red circles indicate when a candidate is being rejected. This is just a sketch of the algorithm; 
circles are not accurately plotted. 

copied cells, and associations between a copied cell and a 
border cell. Both kinds of associations can be extracted 
from the original tessellation. 

After this step, we discard all the old ghost cells and 
replace them with the new ghost cells. When this is done, 
the tessellation algorithm is performed as previously de- 
scribed. This implementation of boundary conditions 
is essentially the same as the method used in AREPO, 
though we require a bit more care because we need to 
retain the tessellation information in the boundary cells. 
As a final note on this formulation, very little of the code 
depends on the number of dimensions; only the tessella- 
tion algorithm itself is significantly affected by D. 



2.7. Mesh Regularity 

For most problems, evolving the mesh according to the 
above prescription can generically lead to cells which are 
long and skinny, and whose mesh generating points are 
very close to their edges. These cells will evolve in an un- 
stable manner, because their faces can move very quickly 
even while their mesh generators are moving slowly, and 
they can also have a very short sound-crossing time. It is 
therefore desirable to steer cells in such a way th at they 




Fig. 6. — The next timestep's boundary cells are identified in 
the x dimension. In this example, we are using periodic boundary 
conditions. The green cells are the interior cells which share a face 
with the boundary. In addition to these cells, we identify their 
neighbors and next-nearest neighbors (blue). All of these cells are 
then copied and moved to become " ghost cells" . 



are using periodic boundary conditions and a cell moves 
off of the computational domain, it is set to be a ghost 
cell and its corresponding ghost cell is set to be a real cell. 
All cells are flagged to be in one of three categories: In- 
side the domain, outside the domain, or a "border cell" , 
meaning it is inside the domain and neighbors a cell out- 
side the domain. 

The next step involves generating a new list of ghost 
cells, first by making copies of all border cells. Additional 
ghost cells get added to this list, by using all neighbors of 
ghost cells which are inside the computational domain. 
This procedure is repeated until the desired number of 
ghost cells is achieved (see Fig. [6|. For our second-order 
methods, we need two layers of ghost cells. 

Next these copies are moved according to the boundary 
conditions, e.g. if the boundaries are reflecting, their 
positions are reflected across the reflecting wall. 

Finally, neighbors are assigned to the copies by using 
the neighbor data from the original tessellation. These 
associations are of two kinds: associations between two 



tend to become more regular. Springel (2010) found 
an effective prescription for this, which is to give the 
mesh generators an additional component to their veloc- 
ity, pointed in the direction of the center of mass. We 
repeat this prescription here: 



Sj-fj di-0.9i)Ri 
di 0.2rjRi 



di/i'qRi) < 0.9 
0.9 < di/^Ri) < 1.1 
1.1 < di/lrjRi) 

(34) 

Ri = \J dVi Jtt is the effective radius of cell i, di is the 
distance between the cell's mesh generating point r*i and 
its center of mass Sj. Cj is the local sound speed, rj and 
X are arbitrary parameters for this prescription, which 
are typically set to r\ = 0.25 and x — 1-0 (the same val- 
ues typically used by AREPO). Note that we we do not 
implement a relativistic velocity addition formula, as w 
need not be interpreted as a physical velocity. 
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TABLE 1 
One-Dimensional Tests 



Test 


Variable 








Nonrclativistic 




x < .5 


x > .5 




Shock Tube 


P 


1 


.25 




N = 100 


V 










r = 1.4 


I' 


1 


.1795 




t = 3.0 










Nonrclativistic 




x < .1 


.1 < x < .9 


x > .9 


Interacting Shocks 


P 


1 


1 


1 


N = 400 


V 











T = 1.4 


p 


1000 


.01 


100 


t = .038 










Easy Relativistic 




x < .5 


x > .5 




Shock Tube 


p 


1 


1 






vx 










N = 400 


vy 





.99 




r = 5/3 


1' 


1000 


.01 




t = 0.4 










Hard Relativistic 




x < .5 


x > .5 




Shock Tube 


p 


1 


1 






vx 










N = 100 


vy 


.9 


.9 




r = 5/3 


P 


1000 


.01 




t = 0.6 











3. TEST PROBLEMS 

As this is a new method for solving relativistic hydro- 
dynamics, we present a large number of test problems 
in one and two dimensions. In addition to the relativis- 
tic tests, we have included some nonrelativistic ones, to 
compare our code with AREPO. 

3.1. One- Dimensional Tests 

All ID tests involving piecewise constant states are 
summarized in table [T] Our first two te sts are identi- 
cal to tests performed by Springel ( 2010[ ). The first is a 
si mple shock tube, a te s t whic h has also been performed 
by Hernquist & Katz (1989). To demonstrate the im- 
portance of the HLLC Riemann solver in our scheme, we 
perform four tests, varying whether the mesh is static or 
moving, and varying whether we use HLL or HLLC. Re- 
sults are plotted in Fig. J7J When the cells were moved 
and HLLC was employecL~the contact discontinuity was 
very well approximated. In the other three cases, there 
was far more diffusion of the contact discontinuity. This 
demonstrates that the accuracy in this scheme comes 
from the combination of using a moving mesh and em- 
ploying a multi-state Riemann solver. 

The second nonrelativistic test we perform involves 
t he int eraction of multiple shocks (Woodward & Colella 

1984 1. Fig. [9] shows the multiple-shock problem at 

= .1)38. We compare results with a fixed and moving 



mesh, and using the HLL and HLLC Riemann solvers. 
Again, we see that the combination of using the HLLC 
solver and allowing the cells to move leads to a very high 
accuracy in the solution. 

For relativistic one-dimensional tests, we compare our 
code t o the relativistic adaptive m esh refinement code 
RAM flZhang fc MacFadyen ]|2006[ ) which was tested for 
convergence against a variety of test problems. The first 
two are Riemann problems with transverse velocity, the 
so-called "Easy" and "Hard" shock tube tests. The easy 
shock tube test is plotted in Fig. [8] at time t = 0.4. For 
the easy shock tube, we find nearly first-order conver- 
gence, and smaller errors than RAM (see table [2] for con- 



TABLE 2 

Li ERRORS OF THE DENSITY FOR THE EASY 
RELATIVISTIC SHOCK TUBE AT t = 0.4. TESS IS 
COMPARED WITH RAM *. 



Code 


N 


L\ Error 


Convergence Rate 


TESS 


100 


4.23c-l 






200 


2.57e-l 


0.82 




400 


1.36e-l 


0.82 




800 


7.45e-2 


0.98 




1600 


3.54e-2 


0.91 




3200 


2.48e-2 


0.95 


RAM 


100 


8.48e-l 






200 


4.25e-l 


1.0 




400 


2.41c-l 


0.82 




800 


1.27e-l 


0.92 




1600 


6.43e-2 


0.99 




3200 


3.34e-2 


0.95 



RAM used piecewise parabolic reconstruction 
and a modified Marquina flux 



TABLE 3 

Li ERRORS OF THE DENSITY FOR THE HARD 
RELATIVISTIC SHOCK TUBE. TESS IS COMPARED 
WITH RAM* 



Code 


N 


L\ Error 


Convergence Rate 


TESS 


400 


7.12e-l 






800 


4.54e-l 


0.64 




1600 


2.26e-l 


1.0 




3200 


1.10e-l 


1.0 


RAM 


400 


5.21e-l 






800 


3.63e-l 


0.52 




1600 


2.33e-l 


0.64 




3200 


1.26e-l 


0.89 



RAM used piecewise parabolic reconstruction 
and a modified Marquina flux 

vergence rates). For the hard shock tube, we find that 
like RAM, we need very high resolution to capture the 
solution accurately. However, because we can initially 
place the cells anywhere we want, we can resolve the 
initial discontinuity very well and capture the sol ution 
with as few as 100 conputational zones (see Fig. fill. 



50 of the zones were concentrated uniformly within a re- 
gion Ax = .005 of the discontinuity and the remaining 
50 cells were distributed uniformly through the rest of 
the domain. Using a uniform initial mesh, TESS showed 
first order convergence for this problem (table p5j). 

The last test we perform in one dimension is to demon- 
strate the convergence rate on a smooth problem. We set 
up an isentropic pulse identical to the one used by Zhang 
& MacFadyen I p006l): 



P = Pref(l 



af(x)) 



\x\ < L 
otherwise 



P = Kp l 



(35) 
(36) 
(37) 



Pref, K,L, and a are constants. In our case, p re f = 1.0 
K = 100, L = 0.3, a = 1. T is the adiabatic index, cho- 
sen to be 5/3. To determine the velocity, we use 



i]n(— ' 
2 y l-v' 



ln( 



Vr- i 



VT^i Vr - 1 



(38) 
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Fig. 7. — Rest mass density in the nonrelativistic shock tube at t = 3.0 with 100 cells. Results for a static mesh on the left, and moving 
mesh on the right. Top plots use the HLL Riemann solver, bottom plots use HLLC. The solid line represents the exact solution. The 
contact discontinuity is much better preserved when employing HLLC on the moving mesh. 




TABLE 4 

Li ERRORS OF THE DENSITY FOR THE ISENTROPIC PULSE. 
TESS IS COMPARED WITH TWO VERSIONS OF RAM. 



Fig. 8. — Easy relativistic shock tube at t = 0.4 with 400 cells 
ising a moving mesh. Solid lines represent the exact solution. 



= constant 

where c s is the sound speed. 

Fig. [10] shows the pulse at t = 0.8 on a mesh with 
400 cells. The LI error and convergence rates are shown 
in table [4j In order to make a reasonable comparison 
between the codes, we chose two different methods em- 
ployed by RAM which are second order. Of course, if 
we had chosen to compare with the WENO solvers em- 
ployed by RAM, it would be no contest, as RAM can 
get up to fifth order convergence for smooth flows. For 
this problem, we not only find smaller errors, but also 
slightly faster convergence than RAM (RAM employed 
piecewise linear reconstruction of the primitive variables, 
and a fourth-order Runge-Kutta time integration, lead- 
ing to an overall second-order scheme). In fact, conver- 
gence for TESS was slightly better than second order. 
This could be due to the fact that the method becomes 
"more Lagrangian" as the resolution increases; face ve- 



Code 


N 


Li Error 


Convergence Rate 


TESS 


80 


4.88e-3 






160 


1.78e-3 


1.8 




320 


4.84e-4 


1.9 




640 


1.20e-4 


2.0 




1280 


2.93e-5 


2.1 




2560 


6.97e-6 


2.1 




5120 


1.47e-6 


2.1 


RAM 


80 


1.12e-2 




(U-PLM-RK4) 


160 


3.56e-3 


1.7 




320 


1.02e-3 


1.8 




640 


2.60e-4 


2.0 




1280 


6.49e-5 


2.0 




2560 


1.62e-5 


2.0 




5120 


4.04e-6 


2.0 


RAM 


80 


l.lOc-2 




(U-PPM-RK4) 


160 


2.56e-3 


2.1 




320 


5.74e-4 


2.2 




640 


1.34e-4 


2.1 




1280 


3.10c-5 


2.1 




2560 


7.33e-6 


2.1 




5120 


1.82e-6 


2.1 



locities approach the velocities of contact waves in this 
limit. 

3.2. Convergence in Multiple Dimensions 

On the surface, it is certainly not obvious that this 
second-order convergence extends to multiple dimen- 
sions; it would be quite difficult to prove such a thing 
mathematically. A convergence test in 2D is essential if 
we want to establish confidence in our numerical scheme. 
A straightforward test which satisfies this is to propa- 
gate an isentropic wave diagonally across the mesh. The 
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Fig. 9. — Rest mass density in the double blastwave test at t = .038 with 400 cells. Results for a static mesh on the left, and moving 
mesh on the right. Top plots use the HLL Riemann solver, bottom plots use HLLC. The solid line is meant to represent the exact solution; 
it was calculated numerically at much higher resolution. Moving the mesh improves accuracy, and using HLLC on this mesh brings a vast 
improvement in resolving both shocks and contact discontinuities. 




Fig. 10. — Isentropic pulse at t = 0.8 with 400 cells; a smooth 
problem to test second-order convergence. The solid line was cal- 
culated at much higher resolution. 

mesh will initially be a square grid, but because the flow 
is nonuniform, the mesh will move in a nontrivial way 
and we will be able to check whether second-order con- 
vergence continues to hold during such distortions. The 
result we present will employ a nonrelativistic wave in a 
periodic box (x, y) £ [0, 1] x [0, 1], with density and pres- 
sure given by the same formulas as in the ID case, but 
now with 



f(x,y) = sin 2 (Tr(x + y)). 



(39) 



For this test, we use K — .01, but all other constants 
are the same as in the p revi ous example. The nonrela- 
tivistic limit of equation ( 39 1 is simply 



1 



r,ref 



(40) 



The direction of this velocity is of course diagonal to 
the mesh, 



v(x + y 



(41) 



Condition (40) continues to hold as the wave evolves, 
so we can use this as a wa y of measuring LI error for 
this problem. In Fig. [12] we see the pulse at t = 1.0 
on a 50 x 50 mesh, and the cells have clearly changed 
their shape and size in a non triv ial way. LI error for 
this problem is plotted in Fig. 13 The convergence rate 



for this problem was calculatedto be 2.3, again slightly 
better than second order. 



3.3. Tests in Spherical and Cylindrical Coordinates 

TESS is written in a modular way so that it is straight- 
forward to change the set of conserved quantities, and to 
add source terms. As an example, we have implemented 
the equations in alternative coordinate systems for the 
two special cases of spherical symmetry and axisymme- 
try. Our extension to these coordinate systems takes the 
most naive approach: we simply express the equations in 
an alternate coordinate system, and solve them like any 
other hyperbolic system in conservation-law form. We 
do not wish to endure the complication of curvilinear 
voronoi cells, which is why we specialize to the cases of 
spherical symmetry (1-D) and axisymmetry (2-D). For 
example, in cylindrical coordinates in axisymmetry the 
nonrelativistic conservation laws take the following form: 



U = 



rp 
rpv r 
rpv z 



r 2 pv<p 
\ r(pv 2 /2 + 



\ 



(42) 
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Fig. 11. — Hard relativistic shock tube at t = 0.6 (solid line is the exact solution). Both tests were run with a moving mesh using 100 
zones, but the lower test concentrated 50 of the zones within a small region near the initial discontinuity, momentarily giving an effective 
resolution of roughly 10,000. 



0.0001 
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Fig. 12. — The 2D isentropic wave at t = 1.0 on a 50 X 50 mesh. 
This test demonstrates second-order convergence when the mesh 
is allowed to move and dynamically change its geometry (initially 
we start with a uniform grid). 



F-h 



( rpv ■ h \ 

r(pv r v ■ h + Pn r ) 
r(pv z v ■ h + Ph z ) 
r 2 pv&v ■ h 
\ r(pv 2 /2 + pe + P)v ■ h J 



(43) 



/ 



S = 



\ 

P 







Fig. 13. — LI error in the velocity for the 2D isentropic wave at 
t = 1.0 at various resolutions. The slope of the solid line is 2.3. 



where r is the distance from the z-axis (note that n lives 
in the r-z plane) . The two new developments here are the 
presence of the radial coordinate r in the equations and 
the non-zero source term. We evaluate both of these 
at the center of mass of a cell or face, depending on the 
context. The description in spherical coordinates and the 
extensions to relativistic hydrodynamics are completely 
analogous. 

We use this opportunity to test the method's exten- 
sion to these coordinate systems, as well as the ability 
for TESS to preserve symmetries using an unstructured 
mesh. We set up very simple shock-tube like initial con- 
ditions: 



(44) 



1.0 
0.1 



r < 
r > 



.25 
.25 



(45) 
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Fig. 14. — Density profile for the relativistic cylindrical explosion 
using a 50 X 50 cartesian mesh at t = 0.2. The plot above uses a 
fixed mesh, while the plot below uses a moving mesh. The solid 
line was calculated in one dimension (cylindrical coordinates) at 
high resolution. The cells plotted lie along a single line; blue pluses 
represent cells that lie along the x-axis and purple crosses represent 
cells that lie along the diagonal x = y. 



P = P 

v = 



(46) 
(47) 



We perform a relativistic cylindrical and spherical ex- 
plosion. For the cylindrical explosion, we perform the 
calculation at high resolution in 1-d using cylindrical co- 
ordinates, and in low resolution in 2D using cartesian 
coordinates. The resolution is as low as 50 x 50. In 
Fig. [14] we show the results for a moving and fixed mesh. 
We see that moving the cells continues to improve res- 
olution of the contact discontinuity and in this case we 
also see some improvement in symmetry - the values for 
the density profile are not as scattered in the moving- 
mesh case, they tend to lie along a single curve. For the 
spherical explosion (Fig. 15), we see a similar story; the 
contact discontinuity is well preserved, and the spherical 
symmetry is captured somewhat better with the moving 
code. 

3.4. Test Problems in Magnetohydrody amies 

To further illustrate that TESS in principle extends to 
an arbitrary set of conserved quantities, we demonstrate 
a few test problems in magnetohydrodynamics. From 
one point of view, MHD (whether relativistic or non- 
relativistic) is just another hyperbolic set of equations 
in conservation-law form, where electromagnetic energy 
momentum, stress, and pressure are added to the fluid 
equations, and the magnetic induction also satisfies a 
conservation law (Faraday's Law), which in the the case 

of a perfect conductor (E = — v x B) has the form 




Fig. 15. — Density profile for the relativistic spherical explo- 
sion using a 50 X 50 mesh in cylindr ical coordinates at t = 0.2. 
Conventions are the same as in Figure 1141 



The main ingredient which distinguishes MHD from 
other conservation laws is that the equations also satisfy 
an elliptic constraint, V • B = 0, which must be satisfied 
every timestep. This constraint is analytically guaran- 
teed to vanish, but depending on your numerical scheme, 
the constraint can generically grow exponentially when 
the equations are solved in discrete form. Over the years, 
many techniques have been suggested to address this is- 
sue, the most pop ular methods being variants of con- 
strained transport (Toth [2000), where steps are taken 
to ensure that all fluxes added to the B fields have zero 
divergence, so that the time derivative of div B is guar- 
anteed to remain zero through the entire evolution (in 
other words, only a solenoidal B field is ever added dur- 
ing the time evolution). This method would be quite 
difficult to extend to an unstructured mesh, so it seems 
that we must choose some other way of preventing these 
constraint-violating solutions from growing. 
We employ a hyper bolic divergence-cleaning scheme 



(Dedner et al. 2002) which adds an extra conserved 



quantity i/j, and modifies the equations so that the evo- 
lution equations for the magnetic field become 



d t B + Vj{vjB - BjV) + Vip = 



ftV + c l V • B 



(49) 
(50) 



d t B + dj{v J B - B j v) =0 



(48) 



where c\ and are freely specifiable constants which 
alter the behavior of the divergence-cleaning terms. If 
V ■ B — and ifj — the equations revert to the usual 
MHD equations, but this introduction of ip has the effect 
of altering the time evolution of the constraint so that it 
satisfies a wave equation, specifically the telegraph equa- 
tion. This technique is generally not as desirable as con- 
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strained transport, partly because it doesn't guarantee 
zero divergence, but also because its success is strongly 
dependent on two freely specifiable constants, which have 
no obvious choice for what their values should be. Be- 
cause we don't have a constrained transport scheme 
(though we do speculate that one is possible), and also 
because this method is easy to implement, we currently 
use this divergence-cleaning form of the equations in our 
MHD code. 

We now look at several MHD tests. The first is a one- 
dimensional shock tube test in relativistic MHD, the rel- 



ativistic version of the B rio-Wu shock tube (Brio & Wu 



1988 
'as tol 



van Putten 1993). The left and right states are 



ows: 



PL = 1: 
Pl = 1, 



PR = -125 
Pr = .1 



(51) 
(52) 



The initial velocity is zero everywhere, and the mag- 
netic field is given by the following: 



B x = 0.5, 



B 7 =0 



D 



yL 



1. 



B 



V R 



1 



(53) 
(54) 



The adiabatic index T = 2. The calculation was done 
with 400 zones. The final state is shown at t — 0.4, in 
Fig. |16| It is clear that an advantage is gained in moving 
the cells, both in the sharpness of the contact discontinu- 
ity and in the shock front. Additionally, the value of the 
density in the shock wave appears to be more accurate 
when using a moving mesh. So, it would seem that there 
is great potential in using this Lagrangian approach for 
problems in magnetohydrodynamics. However, we wish 
to avoid reading too much from this result, as this is a 
one-dimensional problem, and most of the challenges in 
magnetohydrodynamics only appear when solving multi- 
dimensional problems. 

We consider now a test problem in nonrelativistic mag- 
netohydrodynamics, a cylindrical explosion in a uniform 
magnetic field. The calculation is done in a box of dimen- 
sions [0,1] x [0,1]. The initial conditions are as follows: 



P 



P = 
B x - 

10.0 
0.1 



1.0 

= 1.0 



r < 0.1 
r > 0.1 



(55) 
(56) 

(57) 



The magnetic pressure at t=0.1 is plotted in Fig. [17 
Also plotted is the numerically calculated divergence of 
the magnetic field. This divergence is normalized in the 
following manner: 



(divB) no 



V • B 
\B\jR 



(58) 



Where R is the characteristic size of the Voronoi cell. 
We see the code does a reasonable job of resolving the 
explosion at low resolution, however the magnetic diver- 
gence is unacceptably large (the normalized magnetic 
divergence grows to be as large as 0.14). This large 
error in div B appears to be particularly problematic 
for this method. Any time the Voronoi mesh changes 
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Fig. 16. — Density for the Brio-Wu relativistic MHD shock tube 
at t = 0.4 with 400 cells. The upper plot uses a fixed mesh, and the 
lower plot uses a moving mesh. The solid red line is the analytic 
solution. Just as in the nonrelativistic hydrodynamical case, the 
contact discontinuity is captured much better on the moving mesh. 



topology, there can be extremely rapid growth in mag- 
netic divergence due to the appearance of new faces in 
the tessellation. This magnetic growth is so rapid that 
it shows up even in simple test problems like the one 
displayed. In fact, the growth rate is such that con- 
ventional divergence-cleaning techniques for suppressing 
the growth of div B seem to be insufficient to resolve 
the problem. We believe that this problem is not in- 
surmountable, but the eventual solution might require 
a novel approach for treating div B. As such, we have 
much work to do in improving the magnetohydrodynam- 
ics part of our code. Another possible idea is to use the 
same basic numerical scheme, but impose the geometry 
of the cells by hand rather than determining them based 
on the motion of mesh points. If such a scheme were de- 
veloped, constrained transport could potentially be pos- 
sible, since the geometry of the cells could be known at 
runtime. This idea could have many other possible ad- 
vantages, so it is a thought worth pursuing, but for TESS 
as it is currently written, we would like a better means 
for damping the growth of magnetic monopoles. 

On the other hand, for some problems, div B may not 
be the most important concern. If the magnetic diver- 
gence does not have time to grow, we may get to the end 
of a calculation before it starts to have an effect on the 
solution. As an example, we show the relativistic MHD 
rotor test ( Del Zanna et al~||2003 ), a challenging test of 
the robustness ot our numerical scheme. The calculation 
is done in a box of size [—.5, .5] x [—.5, .5]. The initial 
conditions are given by the following: 



10.0 : r < 0.1 
1.0 : r > 0.1 



P = 1.0, 



B a 



1.0. 



(59) 
(60) 
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Fig. 17. — Nonrelativistic cylindrical explosion in a uniform mag- 
netic field at t=0.1. Magnetic pressure is plotted in the left panel, 
while normalized div B is plotted in the right panel. Magnetic 
divergence grows to be unfortunately large; the normalized diver- 
gence is of order 10 percent. 



w 




Fig. 18. — Relativistic MHD rotor test on a 400 X 400 mesh. 
Upper left is density (0.32 < p < 5.8), upper right is pressure 
(p < 3.7), lower left is magnetic pressure (.5-B 2 < 2.3), and lower 
right is Lorentz factor (7 < 1.5). TESS passes this test without 
requiring any fail-safe backup procedures. 

The disc r < 0.1 is initially rigidly rotating with angu- 
lar frequency to = 9.95, so that the the edge of the disc 
has lorentz factor 7 s» 10. The velocity field is zero out- 
side the disc. The calculation was done using a 400 x 400 
mesh, and the results are shown in Fig. at time t = 0.4. 
This is ordinarily a very challenging RMHD test prob- 
lem, but we required no special fail-safe procedures to 
solve it. Actually, for this particular problem, we are at 
a slight disadvantage, since the mesh does not necessarily 
give us adaptive resolution exactly where we want it; the 
rarefaction is under-resolved (we start with an initially 
uniform mesh, but the central cells expand by roughly 
a factor of five in length). However, it seems that, at 
least for our method and for this problem, div B does 
not interfere with our results. 

3.5. Kelvin Helmholtz Instability 

We find that TESS is very good at solving problems 
involving fluid instabilities, particularly in cases where 
a contact discontinuity is being perturbed. The Kelvin- 
Helmholtz instability is one such example, which occurs 
at the interface between shearing fluids. We look at the 
nonrelativistic case first, to compare with AREPO. 



We use a periodic box, [0, 1] x [0, 1], with a stripe in 
the middle half of the box 1/4 < y < 3/4. The pressure 
is uniform throughout the domain, and the density and 
x-velocity are uniform in each region: 



P = 



P = 2.5 

2 : .25 < y < .75 

1 : otherwise 

.5 : .25 < y < .75 
— .5 : otherwise 



(61) 
(62) 

(63) 



This sets up the initial conditions for the instability, 
and to excite a particular mode we add a small pertur- 
bation to the y-component of the velocity: 



v y = w sin(4irx)f(y) 



f(y) = exp(- 



(y -25) 2 

2a 2 



exp(— 



(y--75) ; 

2a 2 



)■ 



(64) 
(65) 



We choose wq = .1 and a = .05/y/2. We use an adia- 
batic index of T = 5 / 3. These initial c onditions are iden- 
tical to those used by Springel] ( |2010[ ). Our visualization 
is different, so it is difficult to directly compare results, 
but it appears that the growth rate of the instability is 
the same for a 50x50 lattice (see Fig. 



19). 



Tests of Galilean invariance, on the other hand, indi- 
cate that our numerics do depend on the moving refer- 
ence frame in which the calculation is performed. When 
adding an overall velocity to the initial conditions, we 
find that the results are significantly changed if the boost 
velocity is large compared to the flow velocity. This can 
be traced back to the fact that the HLLC solver is not 
invariant under Galilean boosts. It would be interest- 
ing to see how differently the code would perform using 
an exact Riemann solver, though this is not a very high 
priority, because TESS is ultimately designed to solve 
relativistic problems. A more optimistic way to inter- 
pret this result is to say that Galilean invariance is not 
a necessary condition for accurately capturing contact 
discontinuities. 

As a relativistic test, we use TESS to calculate the 
growth rate of the Kelvin-Hclmholtz instability in the 
relativistic case. For this purpose, we change the initial 
conditions slightly: the initial density is uniform (p = 1) 
everywhere, the pressure is now P = 1 and the adiabatic 
index T = 1.4. The initial perturbation is given to the 
y component of the four- velocity, u y , and the x compo- 
nent u x is varied. To measure the development of the 
instability, we add an additional passive scalar quantity, 
X, which is advected with the fluid velocity, for which 
we choose X = outside the stripe and X = 1 inside 
the stripe. We measure the growth of the instability by 
tracking the position of the contact discontinuity in X. 
Note that this is much more straightforward to do with 
TESS than it would be for an Eulerian code. In fact, we 
can track the development of the instability even while 
the perturbation to the contact is smaller than the size 
of a Voronoi cell. 

We compare the linear growth rate against the analytic 
solution. The growth rate h as been calculated analyti- 
cally by Bodo et al. (2004). It is the imaginary part 
of (pkc a where k is the wavevector corresponding to the 
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Fig. 20. — Linear growth rates of the relativistic Kclvin- 
Helmholtz instability, calculated using TESS at resolutions of 64, 
128, 256, and 512 (squared). The solid line is the analytical solu- 
tion. 



wavelength of the perturbation (fc = An in this case), c s is 
the sound speed, and is given by the following formula: 



4> 2 _l-v 2 M 2 + l-v 2 - v / AM 2 (l-v 2 ) + {l + v 2 ) 2 
M 2 ~ l-c 2 ~ M 2 + 2v 2 

(66) 

Here, v is the flow speed, and M. is the relativistic Mach 
number, defined: 

M = — (67) 

where c s is the relativistic sound speed, and "f s is the 
correspon ding Lorentz factor. 

In Fig. [20] we plot the numerically calculated growth 
rates against this analytical formula. At 64 x 64 resolu- 
tion, the growth rate was not very accurately calculated, 
but for all other tests we seem to get very accurate re- 
sults, considering the low resolution. 

Notice that in contrast with the nonrelativistic exam- 
ple, it does not make sense to even ask about the Lorentz 
invariance of the code for this problem, since periodic 
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Fig. 21. — Relativistic Rayleigh- Taylor instability, at a resolution of 128x 128, with a relativistic Atwood number of A = 0.6, for a moving 
mesh (above) and fixed mesh (below), plotted at regular intervals until t = 7.0. Asymmetry is primarily the result of small perturbations 
initially given to the mesh generators. These same small perturbations were given in the fixed mesh case, so that the only difference between 
the two versions is whether the mesh points are moving. 
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Fig. 22. — Linear growth rates of the relativistic Rayleigh- Taylor 
instability, calculated using TESS at resolutions of 64, 128, 256, 
and 512 (squared). The solid line is the analytical solution. 

boundary conditions are only periodic in one reference 
frame, due to relative simultaneity. In general, even if 
we do not employ periodic boundaries, testing Lorentz 
invariance can be an extremely complicated task, be- 
cause one would have to take into account that different 
mesh points can be sampling the fluid motion at differ- 
ent points in time. In any case, we know that our code 
is not manifestly Lorentz invariant, so we would not ex- 
pect a striking improvement in Lorentz invariance over 
an Eulerian code. 

3.6. Rayleigh Taylor Instability 

A related fluid instability is the Rayleigh- Taylor insta- 
bility, where a fluid of high density is balanced on top 
of a fluid of low density, and the contact discontinuity 
between these two different densities is perturbed by the 
force of gravity. In the relativistic case, we can get a 
gravitational field by transforming to an accelerated ref- 
erence frame. This is not difficult to do, but we can 
also simply solve the relativistic equations in the weak- 



field limit. For the purposes of calculating the linear 
growth rate, it would not matter which we do, since for 
a small enough perturbation, the potential jump between 
the highest and lowest points of the perturbation will be 
small, and hence the weak-field limit applies. We set 
up the following initial conditions in the computational 
domain (x,y) € [0,1] x [—1,1]: 



1.0 

P2 



y<0 
y>0 



P = P e- 9V/{T - 1] + (r - l)p(e- 9v/{r - 1} - 1) (69) 



For the initial velocity perturbation: 



y 

v y — Wocos(2irx)exp(—-^-^), 



(70) 



We choose constants P = 10, 
and the gravitational field is g - 
different grow th rates. 



w - 

: 0.1. 



= 0.03, cr = .05/72, 
P2 is varied to get 



In Fig. 21 



we show the growth of the instability using 
a moving or fixed mesh, at fairly low resolution, 128x128. 
The initial mesh was not perfectly square; the cells were 
randomly offset slightly to improve the regularity of the 
mesh throughout the calculation. This small ( .1%) offset 
is enough to cause asymmetry, even in the fixed-mesh 
case, but we decided to employ the offset in both cases, so 
that the only difference between the two cases is whether 
the mesh points are allowed to move. 

The linear growth rate of the instability is easy to de- 
rive even in the relativistic case, because the only rela- 
tivistic effect for small perturbations is that gravity cou- 
ples to energy instead of mass. Therefore, the growth 
rate for the relativistic case is still 



R = s/gkA, 



(71) 
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Fig. 23. — Relativistic Rayleigh- Taylor instability in a decelerating shock. The initial density profile is given by a power law, with k = 
(left to right) 0, 1, and 2. In spherical symmetry the initial conditions would result in a transient two-shock solution, but in more than one 
dimension, the contact discontinuity is unstable and it is disrupted enough to generate unstable motion all throughout the forward shock. 
The first two calculations have been terminated at t = .95 when the shock has decelerated to 7 ~ 4, the last one at t = .7, while it is still 
highly relativistic, 7 ~ 11. 



Where g is the gravitational field strength, k is the 
wavevector corresponding to the perturbation (k = 2n 
in this case), and A is the relativistic Atwood number, 
defined as 



A 



P2O- + 62) - pi(l + gi) 
p 2 (l + e 2 )+pi(l + ei) 



(72) 



In Fig. [22] we plot the growth rates numerically calcu- 
lated using TESS, and compare with the analytical result 
(72). Even at extremely low resolutions, we capture the 



growth rates quite accurately. 

4. RAYLEIGH TAYLOR INSTABILITY IN A 
DECELERATING SHOCK 

The Rayleigh- Taylor instability provides us with a use- 
ful astrophysical application, in the context of a decel- 
erating relativistic shock. The reason this is relevant is 
that an external gravitational field is equivalent to an 
accelerating reference frame. Hence a low-density gas 
accelerating into a high-density gas or a high-density gas 
decelerating into a low-density gas will exhibit this same 
instability. The latter could potentially occur in a shock- 
wave, as matter is pushing its way through some ambient 



medium (Levinson 20101. If the explosion is spheri- 



cal, and the density profile has a power law dependence, 
p ~ r~ k , then the shock will decelerate for k < 3. If we 
assume a simple point explosion in this de nsity profile, 
we recover th e Blandford-McKee solution (Blandord & 
McKee ||l976 ), or at late enough times, the nonrelativis- 
tic Sedov-Taylor solution. In both cases, the shock front 



and contact discontinuity coincide. However, if there is 
excess mass in the initial explosion, so that the total 
mass is greater than the integral of the power-law den- 
sity profile, then there will be a coasting period until the 
shock overtakes an amount of mass comparable to this 
excess mass. T his is followed by deceleration to a two- 
shock solution (| Nakamura fc Shigeyama |2006 [) , because 
the information that the shock is now decelerating needs 
to be transported back upstream. In this case, there will 
be a contact discontinuity between the two shocks, and 
the density discontinuity across this contact can be quite 
significant. For this situation, the Rayleigh- Taylor insta- 
bility can play an important role in the solution, because 
the contact is quite unstable. A calculation in spheri- 
cal symmetry would be incorrect, even though the initial 
conditions are spherically symmetric. 

We capture the Rayleigh- Taylor instability from a 
spherically symmetric explosion in two dimensions 
using axisymmetric coordinates, for density profiles 
corresponding to k = 0,1,2. We use the domain 
(s, z) e [0, 1] x [—1,1] where r 2 = s 2 + z 2 and set up the 
initial conditions as follows: 



Po(ro/r min 
P = { Po(ro/r) k " 
Po(ro/r) k 



\k 



r <r m 
r < r 
r > r 



(73) 



We set ko = 4 and ro = .1 to start with an accelerating 
shock which typically coasts until around r ~ .2, then de- 
celerates. We use r m in = -001 to ensure that the system 
has a finite amount of mass. The initial pressure profile 
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Fig. 24. 




A close-up of Figure [23 



is given by 



P 



eo/3 
10" 6 



r < r. 
r > r e 



exp 



(74) 



We use a gamma-law equation of state with F = 4/3. 
We choose r exp = .003, and choose eo to get a shock 
which reaches Lorentz factors of 7 ~ 10 before it deceler- 
ates to mildly relativistic flow by t = .95, when the cal- 
culation ends. The values of eo are given in table[5l along 
with a rough estimate of the bulk Lorentz factorfor the 
fluid at its maximum speed, and also the Lorentz factor 
at the end of the calculation. Note that these initial con- 
ditions ar e fundamental ly different from the conditions 
chosen by |Levinson | pOlO I ; his results came from pertur- 
bative analy s is of t he self-similar solution of |Nakamura fc| 
Shigeyama I ( |2006[ ), whereas ours begins with an explo- 
sion. In spherical symmetry, our explosion approaches 
this self-similar solution, but due to the instability, in 
2D our calculation should never find this solution (at 
late enough times, however, it should eventually find the 
Sedov soluti on) . 

In Figures [23| - [25] we show the growth of the instability 
for three power-law density profiles. It appears that the 
instability grows rapidly enough to catch up to the for- 
ward shock. In fact, relativity appears to help us in this 
case, because for large Lorentz factors the shock front 



TABLE 5 
Rayleigh- Taylor explosion 
parameters 



Power Law (k) 


eo 


'Jmax 


"f final 





6 


11 


4 


1 


1 


9 


4 


2 


6 


13 


11 



must remain a short distance from the contact disconti- 
nuity, keeping it within reach until the Rayleigh- Taylor 
fingers can catch up to the shock. The fact that the insta- 
bility overtakes the forward shock appears to be generic 
and does not require special initial conditions, although 
it would of course not occur for a simple point explosion, 
because in this case the shock and contact discontinu- 
ity would already coincide. Using steep power-laws with 
little deceleration does not seem to change this picture; 
the instability can still catch up to the shock even when 
k = 2 and the lorentz factor is 11, as in our third test 
T his result is in qu alitative agreement with the re- 



casc. 



suits of Levinson (2010), whose work was restricted to 
the linear growth rate. To our knowledge, ours is the 
first direct numerical calculation of the instability for a 
relativistic shock. 
It is worth noting another advantage to TESS which 
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is shown most plainly in Fig. |25| If we know in advance 
where to initially place the cells7 we can get high resolu- 
tion exactly where we need it. In this case, we initially 
concentrated most of the cells near the origin, and the re- 
gion of high resolution followed the motion of the shock. 
Nearly all of the zones have been compressed into the 
shock front, giving us very high resolution there. We gen- 
erated this initial mesh in the following way: We started 
with an initial 512 x 1024 uniform mesh in the domain 
(x, y) £ [0, 1] x [—1, 1] (though we staggered mesh points 
so that our initial mesh was not exactly square) . We then 
changed the location of the mesh points via the following 
prescription: 

r = max{\x\, \y\) (75) 
x -> xe k{r - 1] (76) 

y -> ye k(r - 1} (77) 

where we chose k = 3.5. Thus our initial resolution was 
e k ~ 33 times higher near the origin than it would be for 
a uniform mesh. It should be noted that this effective res- 
olution changes throughout the time evolution, partially 
due to the shock overtaking more cells, and partially due 
to the increasing amount of volume which is being well 
resolved. 

5. SUMMARY 

TESS is a new relativistic hydrodynamics code based 
on the Voronoi tessellation. It performs particularly well 
on problems with sharp discontinuities, especially con- 
tact discontinuities. TESS has been extensively tested. 
For smooth problems, convergence is slightly better than 
second order. For problems involving a moving contact 
discontinuity, TESS demonstrated a clear advantage over 
the Eulerian codes in that smearing of contact disconti- 
nuities is greatly reduced. It was also demonstrated that 
employing the HLLC solver was necessary to get full ad- 
vantage out of the Lagrangian nature of this method. For 



many nonrelativistic problems, TESS has success com- 
parable to that of AREPO, but unfortunately lacks the 
Galilean invariance property. 

We have applied TESS to an astrophysical example, 
testing the stability of a decelerating spherical shock 
wave. We found that the contact discontinuity behind 
the shock was unstable to the Rayleigh-Taylor instabil- 
ity, so much so that the instability was able to reach the 
forward shock. A detailed study will be presented in a 
future publication. 

The structure of TESS is relatively simple, there is 
no need to perform any explicit rotations or boosts, and 
we therefore have the luxury of using an approximate 
Riemann solver, so long as the solver doesn't inherently 
diffuse contact discontinuities. This is a definite advan- 
tage, because for some hyperbolic systems, calculating 
the exact solution to the Riemann problem is computa- 
tionally intensive if not impossible. Our method can be 
used to solve a wide range of hyperbolic problems, like 
general relativistic hydrodynamics, or magnetohydrody- 
namics. The tessellation algorithm is also quite simple, 
and very robust; it takes little effort to implement, and 
is not bug-prone. 

In speed tests, wc find the tessellation algorithm con- 
sumes roughly half of the code's run-time. We urge care 
in the interpretation of this, however, as very little work 
has been spent in optimization, and moreover this figure 
is highly problem-dependent, as is the code's overall run- 
time. On a 2.67GHz Intel Xeon X5550 processor, TESS 
spends roughly 40 microseconds per zone per timestep, 
which is comparable to RAM's efficiency. We hope to 
improve the efficiency in future developments. 

The code is structured in such a way that an exten- 
sion to three dimensions merely requires writing a three- 
dimensional tessellation algorithm, which would be an 
extension of the two-dimensional algorithm already out- 
lined. Making the code run in parallel does not pose 
any major academic hurdles (apart from load-balancing). 
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Both of these ideas are extensions to the code which we 
are currently pursuing. We are also interested in imple- 
menting more complicated boundary conditions, mesh 
refinement, and a local time step, as these have all been 
implemented in AREPO and are thus a solved problem. 
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APPENDIX 

PRIMITIVE VARIABLE SOLVER FOR RELATIVISTIC MAGNETOHYDRO DYNAMICS 

The conversion from conserved variables back to their primitive counterparts is nontrivial for relativistic magne- 
tohydrodynamics. It is particularly challenging if we wish to use an arbitrary equation of state. Our code employs 
a three-dimensional Newton-Raphson solver to recover the pri mitive variab l es fro m the conserved variables. The 
equations being solved are nearly the same as those given by Noble et al. ( 2006 1 ; we derive them again here for 
completeness. 

We begin with the eight conserved variables, which we know: 



D = 7P , B, 
We also make the following convenient definitions: 



fluid 



T 



On 



EM 



K = 7V 



W = ~f 2 ph = VK+1D + {K+ l)(e + P) 
The electromagnetic energy and momentum is straightforward to calculate assuming E — ~v x B: 

T° E ° M = \(E 2 + B 2 ) = x Bf + B 2 ) = B 2 - 7^(B 2 + (u-Bf)h 2 
t em = (Ex Bf = -((« x B) x B) 1 = B 2 v l - B i (u-B)/i 
So, the four-momentum Q can be written out explicitly: 

1 B 2 + (u ■ B) 2 



Q° = w - P + B 2 - 



Q = (W + B 2 )v- {u-B)B/j 
Following Noble, we can also take the dot product Q ■ B to evaluate u ■ B: 



(Al) 



(A2) 
(A3) 

(A4) 
(A5) 



(A6) 
(A7) 



u-B = jB-Q/W (A8) 
If we substitute this back into the equations for the four-momentum, we arrive at the following two relations: 



Q° = W - P + B 2 - -{B 2 /(K + 1) + (B ■ Q/Wf) 
Q 2 = (W + B 2 ) 2 K/{K + 1) - (2W + B 2 ){B ■ Q) 2 /W 2 



(A9) 
(A10) 



Equations (A9) and (A10), along with the definition for W 1A3J), constitute three equations with the three unknown 
variables, W, K, and the temperature T. We can now easily turn this into a Newton-Raphson rootfinding scheme 
where we search the three-dimensional parameter space of W,K, and T for a solution to the three equations. All of 
this assumes an arbitrary equation of state, 



P = P(p,T) = P( ,T) 
WK + 1 



e = e(p,T)=e( 



D 



(All) 
(A12) 



Note that W,K,Td [0, oo), so that we don't have to worry about what to do if one of our variables is too large (this 
would be the case if, for example, we used velocity as one of our unknowns) . If K becomes negative on a given iteration, 
it is reset to zero to prevent taking a square root of a negative, but W and T are generally not constrained. For an 
adiabatic equation of state, we use the same framework but make the definition 



P{ P ,T)=T, e(p,T)=T/(T-l) 



(A13) 
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so that effectively pressure takes the place of temperature as the third variable being solved for. For our initial guess 
values, we use the values of W,K,and T from the previous timestep. However, if the solver fails to find a root, we try 
again using new guess values. One idea is to try an estimate which is exact in the limit that rest mass or magnetic 
energy dominates. In doing so, we obtain the following estimates for K: 

Q 2 (Q ■ Bf 

K ° = J^ KB = D 2 B 2 ( A14 ) 

To take into account both of these possibilities, we use a guess value which is a weighted average of these two estimates: 

DK D + B 2 K B , k . 

K est = D D + B2 (A15) 



If this estimate fails, we can use a method found by Cerda-Duran et al. ( 2008 ) and try a more safe set of guess values, 
based on the maximum possible values (also assuming the Lorentz factor is never larger than 10 4 ): 

e* = (Q° -D-B 2 /2)/D, p* = D, P*=P{p*,e*) (A16) 

W guess = Q° + P* -B 2 /2, A',,,... H)\ T guess = T(p*,e*) (A17) 

To determine when we have converged on the inverted state, we need some error function which will measure how 
close we are to the true solution. For example, one natural choice is 

AW 

S = \^r\ (A18) 

where AW measures how much W has changed since the last iteration. The strategy then would be to iterate until 
S < TOL for some specified tolerance, typically TOL ~ 10~ u . Unfortunately, it's possible for W to be changing very 
slowly while K and T arc still far from the solution. To prevent this from being a problem, we modify the error function 
as follows: 

. .AW. ,Aif, 2 .AT., , Ai . 

<5 = l ir l + ( ir ) +( ^ } (A19) 

This agrees with the previous definition when we are close to the true solution, but away from the solution, it helps 
to prevent "false positives", where 6 becomes small before the true solution has been found. 
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